import numpy as np
import pandas as pd
from sympy import *

data = pd.read_csv('data1.csv')
time = np.array(data['time'])
ti = np.array([int(t[:2]) * 60 + int(t[3:]) for t in time])
x = np.array(data['x'])
y = np.array(data['y'])
# 影长
x0 = np.array([sqrt(xi * xi + yi * yi) for (xi, yi) in zip(x, y)])
# 太阳方位角
omega = np.array([atan(yi / xi) for (xi, yi) in zip(x, y)])
# 太阳直射点纬度
beta = (23.5 * 28 / 91) * 3.1415 / 180
# 北京经度
c2 = 116
# 对经度进行预处理
h = np.arange(1, 4, 0.5)
c1 = np.array([np.array(
    [c2 - asin(-sin(o) * cos(t) / cos(beta)) for (o, t) in zip(omega, np.array([atan(hi / x0i) for x0i in x0]))]) for hi
    in h])
# 求出c1预测值
c1_ans = np.mean(c1)
print(c1_ans)
# 时角
t = np.array([(t - 4 * (c2 - c1_ans) - 720) * 3.1415 / 720 for t in ti])
cos_time = np.array([cos(ti) for ti in t])

flag = 0.001
# 对h进行纬度范围的搜索(步长为0.2)
for hi in np.arange(1, 4, 0.5):
    theta = np.array([atan(hi / xi) for xi in x0])
    cos_time = np.array([cos(ti) for ti in t])
    for psi in np.arange(-90, 90, 0.1):
        wucha_1 = np.array(
            [cos(o) - ((sin(th) * sin(psi) - sin(beta)) / (cos(th) * cos(psi))) for (o, th) in zip(omega, theta)])
        wucha_2 = np.array([sin(th) - sin(psi) * sin(beta) - cos(psi) * cos(beta) * cos_time for th in theta])
        ans = (np.mean(wucha_1) + np.mean(wucha_2)) / 2
        if abs(ans) <= flag:
            # 以当前杆长求出经度
            c1 = np.array([np.array(
                [c2 - asin(-sin(o) * cos(t) / cos(beta)) for (o, t) in
                 zip(omega, np.array([atan(hi / x0i) for x0i in x0]))])])
            print(np.mean(c1))
            print(psi)
            print(hi)

